function [u] = ex32_SolveSystem(dp_lme,s_near,x_nodes,x_samples,w_samples,parameters,options)
% function [u] =
% ex32_SolveSystem_sol(dp_lme,s_near,x_nodes,x_samples,w_samples,parameters,options)
%
% This function assembled the stiffness matrix K and the right hand side
%     rhs corresponding to the cantilever beam problem explained in
%     Timoshenko's book.
%
% Input:
%    p_lme     : shape functions
%    s_lme     : shape functions gradient
%    s_near   : list of neighbors
%    x_nodes   : node points
%    x_samples : sample point
%    w_samples : gauss weigth for each sample point
%    parameters: L (length), D (diameter), nu (Poisson coefficient), E
%                (Young modulus)
%    options   : lme options
%
% Output:
%    u     : vectorial displacement field
%


% Material parameters
E  = parameters.E;
nu = parameters.nu;

nPts = size(x_nodes,1);
sPts = size(x_samples,1);


%% ------------------------------------------------------------------------
% The right hand side rhs is computed
[rhs ind_Dirichlet] = ex32_RHS(x_nodes,options,parameters);


%% ------------------------------------------------------------------------
%  The stiffness matrix is assembled 
nn=0;
for k=1:sPts
  nn = max(nn, length(s_near{k}));
end
nn = min(nn,nPts);

K = spalloc(2*nPts,2*nPts,2*nn*nPts);
C_stiff=E/(1+nu)/(1-2*nu)*[1-nu,   nu,          0 ;...
	                         nu, 1-nu,          0 ;...
	                          0,    0, (1-2*nu)/2];


for ig=1:sPts
  nact = length(s_near{ig});
  B_ig = zeros(3, 2*nact);
  B_ig(1,1:2:2*nact) = dp_lme{ig}(:,1)';
  B_ig(2,2:2:2*nact) = dp_lme{ig}(:,2)';
  B_ig(3,2:2:2*nact) = dp_lme{ig}(:,1)';
  B_ig(3,1:2:2*nact) = dp_lme{ig}(:,2)';
  K_ig_loc = B_ig'*C_stiff*B_ig;

  %assembly
  active=s_near{ig};
  K(2*active(:)-1,2*active(:)-1) = ...
      K(2*active(:)-1,2*active(:)-1) + ...
      K_ig_loc(1:2:2*nact,1:2:2*nact)*w_samples(ig);
  K(2*active(:),2*active(:)-1) = ...
      K(2*active(:),2*active(:)-1) + ...
      K_ig_loc(2:2:2*nact,1:2:2*nact)*w_samples(ig);
  K(2*active(:)-1,2*active(:)) = ...
      K(2*active(:)-1,2*active(:)) + ...
      K_ig_loc(1:2:2*nact,2:2:2*nact)*w_samples(ig);
  K(2*active(:),2*active(:)) = ...
      K(2*active(:),2*active(:)) + ...
      K_ig_loc(2:2:2*nact,2:2:2*nact)*w_samples(ig);
end



%% Dirichlet BCs are applied

%(x=0,y=0)  ux=0 uy=0
ind = ind_Dirichlet(1);
K(2*ind,:)         = 0;
%K(:,2*ind)         = 0;
rhs(2*ind)         = 0;
K(2*ind,2*ind)     = 1;
K(2*ind-1,:)       = 0;
%K(:,2*ind-1)       = 0;
rhs(2*ind-1)       = 0;
K(2*ind-1,2*ind-1) = 1;

%(x=0,y=D/2)  ux=0
ind = ind_Dirichlet(2);
K(2*ind-1,:)       = 0;
%K(:,2*ind-1)       = 0;
rhs(2*ind-1)       = 0;
K(2*ind-1,2*ind-1) = 1;

save('RHS_lme_h5', 'rhs');

%% ------------------------------------------------------------------------
% The system is solved
u=K\rhs;

